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ABSTRACT 

In this paper, I present the calculation of the third and fourth moments 
of both the distribution function of the large-scale density and the large-scale 
divergence of the velocity field, 9. These calculations are made by the mean of 
perturbative calculations assuming Gaussian initial conditions and are expected 
to be valid in the linear or quasi linear regime. The moments are derived for a 
top-hat window function and for any cosmological parameters O and A. It turns 
out that the dependence with A is always very weak whereas the moments of the 
distribution function of the divergence are strongly dependent on Q. A method to 
measure O using the skewness of this field has already been presented by Bernardeau 
et al. (1993). I show here that the simultaneous measurement of the skewness 
and the kurtosis allows to test the validity of the gravitational instability scenario 
hypothesis. Indeed there is a combination of the first three moments of 9 that is 
almost independent of the cosmological parameters Q and A, 

'04)_3(0 2 } 2 ) (9 2 ) 
- 5-^ « 1.5, 

(the value quoted is valid when the index of the power spectrum at the filtering scale 
is close to -1) so that any cosmic velocity field created by gravitational instabilities 
should verify such a property. 
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1. INTRODUCTION 

At large scales, say for scales larger than 10 Mpc, the amplitudes of 
the rms fluctuation observed in the cosmic fields are still below unity showing that 
these fields have undergone only a moderate evolution. The perturbation theory 
is then well adapted to address many problems associated with these scales. For 
instance the growing rate of the fluctuations can be accurately described by the 
linear approximation (e.g., Efstathiou et al. 1988). This is a crucial result for 
the determination of the initial spectrum. However, the statistical properties of 
the large-scale cosmic fields cannot be reduced to the behavior of the fluctuation 
spectrum, as it would have been the case if the fluctuations remain Gaussian. For 
instance the skewness observed in the density distribution at such large scales (see 
for example Bouchet et al. 1993) is in agreement with the theoretical predictions 
based on calculations using the second order perturbation theory (Peebles 1980, 
Goroff et al. 1986, Juszkiewicz, Bouchet & Colombi 1993a). 

The properties of the cosmic velocity field so far have only been considered 
in relation with the density field. In the linear theory, there is a well known 
relationship (Peebles 1980) between the density and the peculiar velocity field, 
u(x) = v(x) - H r, 

u(x) = -/(0,A)^ / J^*(x)d 3 x' (1) 

where /(O, A) is a function depending on the cosmological parameters and <5(x) 
is the matter overdensity field (5(x) = p(x)/p— 1). According to Peebles (1980), 
the function / can be accurately approximated by /(O) = Q and Lahav et al. 
(1991) recently showed that its A dependence was extremely weak. So in principle a 
comparison of the velocity field with the density field should lead to a determination 
of This analysis has been done using the dipole measurement (that gives the 
absolute velocity of our galaxy with respect to the microwave background) and 
assuming that the local density is given by the galaxy field as given by the 1.2 Jy 
IRAS survey (Strauss & Davis 1988) or by the QDOT survey (Rowan- Robinson et 
al. 1990). These analysis lead to the conclusion that O is close to unity (Kaiser & 
Lahav 1989). However this method raised criticisms since the convergence of the 
integral in (1) may be quite slow (see Peacock 1992). This difficulty has been partly 
overcome by Yahil (1991) who uses a local form of the relation (1), 

0(x) = -^V.u(x) = -/(n,A)*(x) (2) 
-no 

to compare the velocity field as determined by Dekel et al. (1991) with the density 
field. Using the IRAS galaxy sample, he also concludes that O should be close 
to unity. However, Shaya, Tully & Pierce (1992) using optical-based study give a 
slightly lower value for / leading to O £ 0.1. In any case, we must have in mind that 
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the light distribution may be a biased indicator for the actual mass distribution. 
For small fluctuations the existence of such a bias can be parametrized by the 
relationship 

5 galaxy( x ) = ( 3 ) 

As a result the methods previously mentionned only give the value of the ratio 
f(Q)/b. It then would be of great interest to find a method that would separate 
these two quantities. Yahil (1991) proposed a way to do it but he assumed some 
arbitrary form for the nonlinear relationship between the galaxy density field and 
the matter density field that has no reason to be true. The nonlinear corrections 
between the matter density field and the velocity field have been investigated either 
by numerical simulations (Nusser & Dekel 1991) or on theoretical basis (Bernardeau 
1992b). However these corrections do not give a way to separate the measurement 
of from the bias effects (see section 4.2 of this paper). 

One way to overcome the problem of galaxy biasing is simply to use only the 
velocity field as a probe of the large-scale structures. The linear theory (eq. [1] 
or equivalently [2]) is obviously of no help. However some statistical properties of 
the velocity field can be derived from perturbative calculations at higher orders. It 
is noticeable that the divergence of the peculiar velocity field is a scalar field that 
can be analyzed in a similar way than the density field. And as data on large-scale 
velocity field begin to be available (Lynden-Bell et al. 1988, Dekel, Bertschinger 
& Faber 1991) it becomes crucial to explore the various statistical indicators that 
could be of cosmological interests in such a field. Following the work of Dekel 
et al. (1991) it seems possible to build the full 3D velocity field from the mere 
knowledge of the line of sight velocities of a sample of galaxies. The reconstruction 
is based on the assumption that the large-scale velocity field must be only potential 
- non-rotational - so that it derives from a simple scalar field. The divergence of 
the velocity field then contains as much information as the peculiar velocity field 
apart from a possible uniform flow which just corresponds to the invariance of the 
equations of motion under a Galilean transformation. 

In the following I will assume that the large-scale structures formed by 
gravitational instabilities from initial Gaussian fluctuations. There are indications 
that the initial fluctuations were actually Gaussian in the galaxy distribution (see 
further in this introduction). Then quantities such as the skewness, (5 3 ) and (# 3 ), 

2 2 

or the kurtosis, (5 4 ) — 3(5 2 ) and (6> 4 ) — 3(6> 2 ) , of the cosmic fields are then of 
great interest ((.) is an ensemble average over the initial fluctuations). Indeed, 
because of the general properties of Gaussian variables, these quantities are zero for 
the initial fields, and the use of the linear theory would conclude that they remain 
zero. But in fact gravity, because of the nonlinear it ies it contains, induces specific 
non-Gaussian features that, even at very large-scales, will be exhibited in these 
parameters in a very precise way. The first example of such a calculation was given 
by Peebles (1980) for the skewness. He showed that its leading term for small a 
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can be derived from a second order perturbative calculation (and higher orders give 
only negligible corrections). Fry (1984) extended these results and showed that the 
kurtosis can be obtained from third order calculations. These results were extended 
afterwards in several directions. A method was presented by Bernardeau (1992a) to 
get the whole series of the cumulants. The O dependence and the smoothing effects 
have been calculated for the density skewness (Bouchet et al. 1992, Goroff et al. 
1986, Juszkiewicz, Bouchet & Colombi 1993) and a recent paper (Bernardeau et al. 
1993) extends these results to the velocity divergence presenting the dependence 
of its skewness with O and the shape of the power spectrum for various window 
functions. The result reads, 



<* 3 > „ -i 



w 2 



o°- 6 



y-(3 + n) 



(4) 



for a top-hat window function (volume-weighted filtering) and a power law 
spectrum of index n. As can be seen in the relation (4) the skewness has 
a dependence with the variance -proportional to its square- that is specific of 
Gaussian initial conditions. The generic relation for non-Gaussian fluctuations 

would have been (# 3 ) oc (6 2 ) 3/2 . Such a property is also true for the density field 
and for the same reason. It has been successfully checked in the IRAS galaxy 
survey (Bouchet et al. 1993), bringing strong support to the Gaussian fluctuations 
hypothesis. The relation (4) has been proposed to determine O from the observed 
velocity field (Bernardeau et al. 1993). The purpose of this paper is to extend the 
result (4) to non-zero values of A and to derive a similar property for the kurtosis 
for both the density and the velocity fields. In part 2, I present the principle of the 
calculations and in part 3, the results. They are restricted to a top-hat window 
function due to its special geometrical properties presented in the appendix. Part 
4 is devoted to comments and comparisons with numerical simulations. 



2. THE METHOD 

2.1. The initial conditions 

At large scale the divergence of the velocity field, as the density field, has 
fluctuations, the distribution of which can be examined by perturbative calculations. 
The first order solution in <5(x) is directly related to the first order in 0(x) through 
the continuity equation (Eqs. [1-2]). The variance of this field is then, in the linear 
regime, driven by the first order approximation. It reads 

<£(x) 2 ) = /(0) 2 <5(x) 2 > (5) 

when the growing mode only is taken into account. 
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In the following, I assume initial Gaussian fluctuations. The density field, at 
first order, can then be written 



^)(x,t) = D 1 (t)e(x) = D^t) f ^!^ k e ik - x (6) 

where D\(t) if the growing factor (see part 3 for the equation the solution of which 
Di(t) is) and the random variables e^. follow the rules, 

< e k e k') = Mk + k')P(fc); 

< e k 1 --- e k 2p+1 > = °; (7) 

( e ki ■ ■ ■ e k 2p ) 2 p pl ' ^ e k si e k S2 ) ' ' ' ( e k S2p _ 1 e k S2p ) ■ 

permutations (s) 

where the function 8d is the 3D Dirac distribution, and P(k) is the spectrum of the 
fluctuations. The divergence of the velocity field now reads, 

fl«(x,t) = -f(n)Di(t) J j£^^ k - x > (8) 

at the first order in ej^. At this order, the fields (6) and (8) are purely Gaussian and 
their whole statistics is determined by the shape of the power spectrum. Specific 
features induced by gravity appear when terms of greater orders are taken into 
account. 

2.2. The equations of the dynamics 

Throughout the paper the matter field is approximated by a nonrelativistic 
and self-gravitating fluid at zero pressure. I do not exclude the existence of a non- 
zero cosmological constant A. The expansion factor, a(t), is then solution of the 
equation, 

?=-^ + * (9) 
a 3 3 

For convenience I will characterize the (time dependent) cosmological parameters 
by Q and by the Hubble constant, H, and its time derivative. Let me recall that Q 
is defined by 

AivGp=—H 2 . (10) 

(At the present time H — H .) The fluid is described by the overdensity field 5(x, t) 
and by the peculiar velocity field u(x, t) assumed to be non-rotational so that 



V x x u(x,£) = 0. 



(11) 
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The equations of the dynamics are then given by (e.g., Peebles 1980), 

-o(x, t) + - V x • [(1 + 5(x, t))u(x, t)} = 

|-u(x, t) + -u(x, t) + -(u(x, t).V x )u(x, t) = -- V x t) (12) 
ar a a a 

V x V(x,t) = 47rGpa 2 5(x,t), 

where the spatial derivatives are taken with respect to x. The first step is to take the 
divergence of the second equation so that the system now reads, using the definition 
of O, 

-U(x, t) + [l + o(x, t)]0(x, t) + 4r Vx .[5(x, t)u(x, *)] = 
il ail 

2 + j (9(x, t) + l(9(x, *) + o(x, *) + -Lv x -[u(x, t).V x ]u(x, *) = 

(13) 

where a dot is for a time derivative. 

It is convenient to introduce the Fourier transforms of these fields, 5(k, t) and 
6>(k, £), defined by 

/d 3 k 
^-y3^o(k,t)exp(ix.k), 

0(M) = J -A^#(k,t)exp(ix.k). 
Due to the absence of rotational part in the velocity field, u(x) simply reads 

U(X ' t) = / (2^)3/2 " fc 1 2 k ^ k ^) eX P( iX - k )- ( 15 ) 

Using the relations (14-15) I can now write the equation of the dynamics in the 
k-space: 

^o(k, t) + 0(k, t) + J V(k\ k - k') S(k - k', t) 6(k\ t) = 

2 + 0(k, t) + 1 d(k, t) + ? Oo(k, f) = 

/ (2^7* ^ (k '' k " k ' } + P(k " k '' k ' } " 2S(k '' k " k ' } ] ^ (k " k '' t] 6{k '' t} ' 

(16) 
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where 



V(k, k') = 1 + 



k.k' 



(17) 



Q(k,k') = l 




(18) 



k 2 k' 2 ' 



These functions will be extensively used in the following. The system (16) is the 
basic system of equations used for the perturbative calculations. 

2.3. The filtering of the Gelds 

In practice we have to consider the density and the velocity field after they 
have been filtered by a particular window function. Indeed perturbative calculations 
are only valid at (very) large scale, so that the strong small-scale fluctuations should 
be smoothed out. I make the assumption that the existence of nonlinear fluctuations 
at small scales will not change the result. This assumption has been verified for 
instance for the time dependence of the variance at large scale in various numerical 
simulations (see for instance Efstathiou et al. 1988). There are no reasons to think 
that this is not true for the purpose of these calculations and the comparison with 
numerical simulations presented in part 3 confirm the validity of this hypothesis. 

For technical reason I only consider the case of the top-hat window function. 
However, nothing in practice prevents to use any other window function. For 
instance Goroff et al. (1986) present results for the density field using a Gaussian 
window function but in such a case the smoothing corrections have been to be 
calculated numerically for each power spectrum. For a top-hat window function all 
the results can be given analytically as shown in the next part. 

As it has been noticed by Bernardeau et al. (1993) it is equivalent to apply 
the filter to the velocity field or to its divergence: both are linear operations that 
commute each other. Then I focus my interest on the filtered density and divergence 
of the velocity field. I define S(Rq) and 9(Rq) by 



for a scale Rq. Note that, for the velocity field, it corresponds to a volume weighted 
filtering: the velocity is equally weighted in each point of space, regardless of the 




(19) 



(20) 



with 



W XH (x) = 1 if |x| <R ; 
M^jj(x) = otherwise; 



(21) 
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local density. In practice, the velocity is known only at the positions of the galaxies 
(or the matter points for a numerical simulation) so that it is necessary to smooth 
the field in two steps: the local field has to be defined in each point by a small scale 
(mass-weighted) filtering and then filtered with a volume weighted procedure. This 
is what have been done by Juszkiewicz et al. (1993b). In (21) the filter is applied 
at the origin but the statistics of the filtered fields obviously does not depend on 
the position at which the filter is applied. I also make the hypothesis that a volume 
average is equivalent to an ensemble average so that the moments of distribution 
functions of 5 and 9 are in principle measurable in a large enough sample. 

The quantities (21) can be expressed with the Fourier transforms of 
respectively the density field and the divergence field. They read, 

6{Ro) = / (2^^TH( k ^o)5(k), (22) 
d ( R o) = / ^^ W TH( k ^o) W, (23) 
where Wpjj(fc-Ro) is the Fourier transform of the top-hat window function, 



3 

VF TH (/ci? ) = „ , 3 (sm(kR ) - kR cos(kR )) . 



The variance of these smoothed fields is then given by, 

a2{Ro) = I j£w~ 2 WW WTR{k Ro) < 5(k)5(k/) > w ™ {k> Ro) 
= Dl{t) I J^ w ™ {k ro) m ' (24) 

al(Ro) = / 2 (0) Dl(t) J ^W^{k R ) P(k). (25) 

These results are valid at large scale where, by definition, the linear approximation 
can be applied to the determination of the variance of the density field. The validity 
domain of the results given in the following is in principle the same as the validity 
domain of (24-25). 



2.4. Calculation of the skewness 

The principle of this calculation has been presented in previous papers (Goroff 
et al. 1986, Bouchet et al. 1992 for the density field, Bernardeau et al. 1993 for the 
divergence of the velocity field). The aim is to calculate the leading order of the 
third moment of the distribution functions of S(Rq) or of 9(Rq). The calculation is 
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based on an expansion of these quantities with respect to the random variables 
that describe the initial density and velocity fluctuations. So I write 

S(Rq) = (Ro) + 6^ ( Ro ) + 5 (3) (i? ) + . . . (26) 

where J^ 1 ^ (-Ro) is proportional to e^, 5^ 2 \Ro) is quadratic in and 5^ 3 \Ro) is 
cubic in ejj.... A similar expansion can be written for 9(Rq), 

6(R ) = e^\R ) + e^\R ) + e^(Ro) + ... 



(27) 



3 3 

Then the ensemble averages of [6{R )}° and [6(Ro)]° read 



> 



(5(R ) 3 ) = < [S {1) (Ro) +5 {2) (Ro) +5 {3) (Ro) + • 

= < [s™ (Ro)] 3 > + 3 < [s™ (Ro)] 2 d ( V (Ro) > + .... 
(9(R ) 3 ) = < [o {1) (R ) + (2) (#o) + e^(Ro) + . . .] 3 > 

= < \e^(R )] 3 > + 3 < \eM{Roj\ 2 e^(Ro)) + ... 



(28) 



(29) 



The first terms of these expansions are zero due to the general properties of 
Gaussian variables. The leading terms will be respectively 3( [f^^ 1 ^ (Ro)] ^ 2 KRo)) 

and 3( [^-^(R^)] 9^ 2 \R )) that involve the product of four Gaussian variables. 
Then the ratios, 



5 3 = 



(S(Ro) 3 ) 
(S(Ro) 2 ) 2 
(0(Ro) 3 ) 



2 ' 



(30) 



(31) 



are finite at large scale. The purpose of the next part is to calculate these limits £3 
and S30 as a function of the cosmological parameters O and A and the shape of the 
power spectrum. 

2.5. Calculation of the kurtosis 

The skewness is a sort of measure of the asymmetry of the distribution 
function. The kurtosis measures the flatness of the distribution function compared 
to what would be expected from a Gaussian distribution. The leading order of this 
term involves both the second order and third order of the dynamics. It reads 

(d 4 (R )) - 3(S 2 (Ro)) 2 ~ 6 < [s^(Ro)] * [s {2) (Ro)] ' > + 4 < [<5« (Ro)] * d^(R )) 

-6< [^ORo)] 2 ) < [5^(Ro)] 2 )-12([5^(Ro)] 2 ) (6™ (R )5^ (R )) 

(32) 
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for the density field and 

(e\R ))-S(e 2 (R )) 2 ^ 6< [^(^o)! 2 \eW(R )} 2 ) + 4( [^(ifo)! V 3 >(i?o)) 



-6([e^(R )\ )([e {2 \Ro)\ )-12([e^(Ro)\ > (e^(R )e^(R ) 

(33) 

for the divergence of the velocity field. In the calculation of these differences, we 
can notice that the terms that are subtracted actually appear in the expression 
of the fourth moment. They contain the non-connected part of the moment (so 
called since they factorize in products of moments) and will exactly vanish when 
the fourth order cumulant is calculated. For instance the fourth cumulant of the 
density formally reads 

<5 4 > c ^6< [^(tfo)] 2 [5 (2) (^o)] 2 ) c + 4< [5^(Ro)] 3 S^(Ro)) c , 

where the subscript c means that only the connected terms are retained. We call 
then £4 and S40 the large-scale limits of the ratios, 



S 4 



S49 



(5 4 )- 3(5 2 Y 

(e 4 )-3(e 2 y 



(34) 
(35) 



The next part is devoted to the derivation of the expressions of S3, S30, £4 
and S^q. 



3. THE CALCULATIONS 

3.1. The dynamics of the spherical collapse 

To calculate the previous limits, the dynamics of the spherical collapse is 
quite helpful. In general the equation describing the evolution of the size of spherical 
symmetric perturbation reads, 

k= _GM£ R ^ (36) 

where R is the size of the perturbation and M(< R) is the constant mass within 
the radius R. The dot stands for a time derivative. This equation can be written 
in term of the overdensity 5 defined by, 5 = (R/Rq) — 1, where Rq is the initial 
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(comoving) radius of the perturbation. The equation (36), is equivalent to (using 
Eq. [10]) 

1 2 1 2 HS 1 6 4 5 2 1 2 

"2 °ii^"3(n^ _ 3(i^^ + 9(n^""2 (37) 

The latter equation can be linearized in 5. Let me write this linear term as 

5^(t)=D 1 (t)5 t 

where Si gives the strength of the initial fluctuation and is supposed to be small. 
The function D\ (t) is the time dependence of the growing mode and it is the solution 
of the equation, 

D 1 +2HD 1 - -H 2 flD 1 = 0, (38) 
2 

which is a growing function of time. The function Di is proportional to the 
expansion factor when t is small or for an Einstein-de Sitter universe, but this 
is not true in general. The second order in 5 (when Si is assumed to be the small 
parameter) reads 

5^(t)=D 2 (t) 5 l 
where the function D 2 (t) is solution of the equation, 

D 2 + 2HD 2 - -H 2 fi L> 2 = 3# 2 O D\ + -Di, (39a) 
2 3 



and verifies 

21 



34 

D 2(t) ~ 7^D 2 (t) when t -> 0. (396) 



Similarly the third order reads 



with 

D 3 + 2HD 3 --H 2 ttD 3 = 9H 2 QD 1 D 2 + 8D1 D 2 ~ 8£>i Di (40a) 
2 

and 

682 o 

D 3 (t) ~ — when t -> 0. (406) 

toy 

In the following I will use the functions Di , D 2 and D 3 and their time derivatives 
to build the second and third order of the density and velocity fields. There are two 
other functions that are convenient to consider. They are given by the evolution 



of the divergence of the velocity field, given by 9(t) = 3R/HR — 3 all along the 
collapse. This function can be expanded with respect to Si in the form, 



0(f) = -/(ft, A) E 1 (t)S l + E 2 (t) 6 f + E 3 (t)^ + ... , (41) 



Where the factor f(Q, A) is given by definition by 



/(0,A) 



a d£>i 



(42) 



D x da ' 



From the equation of the dynamics we get 



E 1 (t) =D 1 (t); 




3D 2 D 1 + 2D\] . 



(43) 



(44) 



The exact analytic form of these functions can be obtained from the solution 
of the spherical collapse (e.g., Peebles 1980), but they are complicated and 
it is unnecessary to give them explicitly. For a non-zero A they have to be 
integrated numerically from the equations (39-40). The time variation of the ratios 
D 2 (t)/D 1 (t) : -D 3 (t)/.Di(t) 3 ,... appears simply as a variation of these ratios with 
the cosmological parameters (O, A) as they cover a given trajectory in the (O, A) 
plan. Here, I present the results of numerical integrations only in the case of zero 
curvature, Q + A/3H 2 = 1 and the analytical results for A = 0. The ratios are given 
as a function of O in Figs. 1-2. 

3.2. Expressions of8^(R ) and 0^(R ) 

The fields at the first order are given by the linear solution (Eqs. [6, 8]). As 
a result we have 




(45) 



and 



6^(R ) = -f(^A) J j^^e k W TR (k R ) D^t) 



(46) 



3.3. Expressions of5^(R ) and {2 \R ) 
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We have to write the system (18) at the second order in e k . It leads to the 



system, 



±tf a )(M) + * (a) (M) = «^A J ^^(k',k-k')e k _ k ,e k 
' 2 + §2] ° {2) ^ *) + ^ ^ ( k ' *) + \ ^ (k, t) = 



da 



a- 



3i / 



d°k 



(2tt) 3 / 2 



[P(k', k - k') + P(k - k', k') - 2Q(k', k - k')] e k _ k , e k , 



(47) 



It is convenient to define the notations Vij and Qij by 



Vij — P(kj, kj), Qij — Q(kj, kj). 



In the following the indices will run between 1 and 3. As a result of Eq. (47) 
we obtain 



d 3 ki d 3 k 2 

(27r) 3 / 2 (27j3/ : 



W TH (|ki + k 2 | i? ) 



+ 2 1)2 2l > 2 



(48) 



and 



3 3 

^ (2)( ^ o) = " / {L)W ^W ^TH (Iki + k 2 | Ro) 



xf(Q,A)[D 2 1 



+ ^ E 2 Q.1,2 



(49) 



These results have already been obtained by Bouchet et al. (1992) for the density 
field and by Bernardeau et al. (1993) for the velocity field in case of A = 0. 



3.4. Expressions of5^(R ) and # (3) (#o) 

The calculation of these functions is based on the same principle as for the 



-14- 



second order. The system (16) has now to be written at the third order in e^, 



1 • 



H 



I 



o (3) (M) + # (3) (k,t) 



d 3 k' 

(2tt) 3 /' 



; V{k', k - k') \SW (k - k') 6^ (k') + 5™ (k - k') 6™ (k') 



24b (3) (M) + ^ (3) (M) + !o^(M) 



H 2 

d 3 k' 

(2tt) 3 / 5 



H v ' ' 2 

[P(k',k-k')+P(k-k',k')-2Q(k',k-k')] 6»( 1) (k-k')6»( 2 )(k / ) 

(50) 



I have to define an extra function that will be useful for the intermediate 
calculations, F 3 (t), defined by the equation, 



f 3 + 2hf 3 - hi 2 n f 3 = ^h 2 n D\ 



(51) 



and which behaves like 9 Df(t) /10 when t is small. This function will not enter in 
the expressions of the kurtosis so this equation has not to be solved. Then we can 
search a solution of the system (50) having the form, 



^(k) = (27r) 3 / 2 



d 3 k! d 3 k 2 d 3 k. 



^k/k/ks Mk-k!-k 2 -k 3 ) 



(2tt) 3 / 2 (2tt) 3 / 2 (2tt) 3 / 2 ' 
x [Df Tl^ku k 2 , k 3 ) + £>i D 2 K 2 (k u k 2 , k 3 ) 
+D 3 K 3 {k u k 2 , k 3 ) + F 3 K 4 (k u k 2 , k 3 )] 



(52) 



6> (3) (k) = (2tt) 3/2 



x 



d 3 kx d 3 k 2 d 3 k 3 
(2tt) 3 / 2 (2tt) 3 / 2 (2tt) 3 / 2 
dD 1 ^ 2 n ly _ , , , d£>i 

da 



e k x e k^k z M k - ki - k 2 - k 3 ) 



a—r- F)\ 5i(ki, k 2 , k 3 ) + a^ 1 — - D 2 <S 2 (ki, k 2 , k 3 ) 



da 

+a^-^D 1 S 3 (k 1 , k 2 , k 3 ) + a^^- <S 4 (ki, k 2 , k 3 ) 
da da 

dF 3 



(53) 



+a— — <S 5 (k!,k 2 ,k 3 ) 
da 

where the functions TZi . . . S$ are homegeneous functions of the wave vectors 
(ki, k 2 , k 3 ) and have to satisfy the system (50). From the properties of the functions 
Di, D 2 , D 3 and F 3 (Eqs. [38, 39, 40, 51]), one can find a natural solution of this 
system. Its writing is a bit complicated so that I introduce new simplified notations, 



Vij,k = V{ki + kj, k fc ), Vijk = V(ki, kj + k k ). 



and similar definitions for 2r/,fc and Qijk- Note that Qij,k = Qk,ij but in general 
^Pij,k 7^ Pfc,£? • The homogeneous functions IZi, . . . ,S$ then read 



S3 
S 5 



^2 

P-3 = o 23,12 2l,2! 
*$2 = T (~ Pl2,3 + 3 63,12) Ql,2! 



^ [^3,12 + Pl2,3 — 3Q342] Ql,2! 



2 / 1 1 

g 23,12 Pi, 2 — ( 2^3,12 + g 23,12 ) Ql,2! 



11 \ /5 3 15 

2^3,12 — 2^12,3 + 23,12 j Pi, 2 + ( 7^12,3 + 2^3,12 2~^ 3 > 12 ) ^!' 2 ' 



- (— ^3,12 + 3 23,12) 2l,2! 

-K 3 ; 
-K4. 



So we now have the expression of 8^'(Ro) and 6^ 3 \Ro) (using eventually the 
functions E2 and E% for the latter). These expressions seem rather complicated 
but further simplications are unnecessary. The functions P and 2 defined in (17- 
18) are well adapted for a convolution with a top-hat window function as shown in 
the appendix. 

3.5. Expressions of the skewness and the kurtosis 

The expressions (28) and (32) are then calculated using the general property 
(7) of Gaussian variables. It leads to the expressions for £3 and £4, 



^ =4 / $^§T* P{kl) P{h) Wl w ' 2 WlM * <64) 
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and 



.S^ (.Ro) = 



6 



d 3 ki d 3 k 2 d 3 k 3 



P(fci) P(k 2 ) P(k 3 ) W x W 12 W 23 W 3 



(2tt) 3 (2tt) 3 (2tt) 3 

X [^1,2^2,3 +^1,2^3,2 +^2,1^2,3 + 1^2,1^3,2 

d 3 ki d 3 k 2 d 3 k 3 



+ 4 



(2tt)3 (2tt)3 (2tt) 3 
x [D 3 n x + D 1 D 2 U 2 +D 3 TZ 3 + F 3 1Z 4 ] 



P(h) P(k 2 ) P(k 3 ) W 123 W x W 2 W 3 



(55) 



where 



W t = W TR {ki Rq), Wij = W TR (|k, + kj\ Rq) , 

W ijk = W TR (|ki + kj + k fc | R ) . 

Similar expressions can be obtained for the velocity field. The previous calculations 
lead to technical results that fortunately reduce to very simple expressions. Indeed 
the values of the integrals appearing in these expressions are all given at the end of 
the appendix and are very simple. Careful summations of these terms (that must 
take into account all symmetry factors) lead to the following results, 



S 3 (R ) 
S 3 e(Ro) 

S^Rq) 
£40(-RC)) 

where 



D 2 
■1 



3 7^2 +7i 



/(n,A) 

D 3 

4 4 +12 

1 

/ 2 (0,A) 



3 T^2 +7ie 



E/2 

D\ + 
E 3 



12- 



D2 
Dl 

E 2 



U Tt> " 2 ) 7i + ^7i + ^72; 



Df Df 



12 



Do 



7i = 



lie = 



dlog[a 2 (P )] 
d log R 

dlQg[<7g(flo)] 
d log R 



Dl ' D\ 

d 2 \og[a 2 {R 



lie 



;lle 



(56) 
(57) 
(58) 
\l20m 



12 = 



129 = 



dlog R 
d 2 log[a 2 (i? ) 
dlog 2 P 



In all cases only the leading term has been given. The expected corrections, 
due to higher order terms, contain an extra factor <j 2 (Rq) (or cr 2 [.Ro]) so they are 
expected to be negligible at large scale. 
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In Figs. (1-2) it is shown that the O and A dependences of D 2 /Dj : Ds/Df, 
E<ilD\ and E^/Df are extremely weak and, for the values of Q and A of interest, 
can be safely approximated by their value when = 1, A = 0. So we end up with 



S 3 (#o) = y + 7i; 



^(-Ro) = 
S4e(Ro) 



-1 



60712 

1323 
1 



26 

y +1ie 
62 



3-71 + g-71 + 372; 



+ 

"12088 
441 



+ 



338 



2]~7i0 + ^fiie + g720 



(60) 
(61) 
(62) 
(63) 



One can notice that the effects of smoothing introduce a dependence of these 
coefficients with the scale through the successive logarithmic derivatives of the 
variance. If the smoothing were not taken into account one would have obtained 
S3 = 34/7 and S4 = 60712/1323 (e.g., Bernardeau 1992) and it corresponds to a 
scale independent cr, not necessarily to a small smoothing radius limit. Moreover 
34/7 does not appear as an upper bound for the skewness of the density field. The 
smoothing actually mixes scales and for a variance decreasing with scale it lowers 
both the skewness and the kurtosis. However for particular initial power spectra 
(such the ones given by Hot Dark Matter models) the skewness may take values 
greater than the one obtained without smoothing. Such a trend is indeed observed 
in numerical simulations (Bouchet & Hernquist, 1992). 



3.6. Power law spectrum and scale dependence 

The scale dependences that appear in the expressions (60-63) are only due 
to a possible change of the shape of the power spectrum with scale. A self similar 
power spectrum, i.e. , a power law of k, 

P(k) oc k n , (64) 

leads to values of S3, S30, S4 and S40 that are independent of scale (as long as 
the quasilinear regime is valid). The previous results are then very simple since 
72 = 726* = 0, and read 



34 

S 3 = y-(n + 3); 



S39 = 

S4 = 
S49 = 



-1 

60712 

1323 
1 



26 

y 



(n + 3) 



-y(™ + 3) + 3> + 3) 2 ; 

12088 338, , 7, n2 
liT- 2T(»+3)+3(n + 3) a 



(65) 
(66) 
(67) 
(68) 



In such a case the results for the skewness have already been given in previous 
papers (Juszkiewicz et al. 1993, Bernardeau et al. 1993). 

These functions are given in Figs. 3-5 as a function of Q when A = and 
when + A/ 3H 2 = 1 and for different values of the index of power law spectra. In 
Fig. 6, I present these functions as a function of the power law index n. In practice 
it is very useful to approximate the actual scale dependence of the variance of the 
density field by a power law, 

a 2 (R ) oc R- (3+n) (69) 

where 

dloga 2 OR ) 

n = 

d log #0 

With such an approximation one gets the relationships (65-68). For the spectra of 
interest it turns out to be a good approximation (see Figs. 7-8). In these figures I 
give the expected behavior of the skewness and the kurtosis as a function of scale 
using eqs. (60-63) (thick lines). I also give the kurtosis using the approximation (67) 
and the relations (67, 68) (thin lines). The decrease of the skewness and the kurtosis 
is due to the variation of the index n with scale (upper panel). This variation is 
more important for the CDM spectrum (Fig. 7) than for the spectrum proposed 
by Peacock (1991) (Fig. 8) to fit the observation in the galaxy distribution. It 
corresponds to the well known discrepancy between the shape of the CDM spectrum 
and the one derived from the actual galaxy clustering properties. 

In table 1, I summarize the results of these perturbative calculations for 
the skewness and the kurtosis for an Einstein-de Sitter Universe. The parameters 
obtained in this paper for the density field are somehow lower than the ones obtained 
by Goroff et al. (1986) but this is due to the fact that they used a Gaussian window 
function instead of a top-hat window function. This effect is particularly important 
when the dependence of a with the scale is strong. However, the dependence of the 
results with the cosmological parameters, O and A, which has not been considered 
by Goroff et al. , is roughly independent of the shape of the window function. 

3.7. Validity domain and comparison with numerical simulations 

The results that are given throughout the paper are thought to be valid in 
the linear or quasi-linear regime. Numerical simulations have been done to check 
the evolution of the variance compared to what is given by perturbative calculations 
and the agreement with the predictions based on the linear approximation hold for 
a 2 up to 0.8 at a 15% precision level (e.g., Efstathiou et al. 1988 or Hamilton et 
al. 1991). In principle the results given in this paper are valid in the same domain. 
Indeed, results of numerical simulations show that the agreement for S3 and S4 
(Bouchet & Hernquist 1992) are extremely good even for values of a close to unity. 



(70) 
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Table 1: Predictions of the perturbation theory for S3, S4, S30 and S40 
vs. spectral slope for an Einstein-de Sitter Universe. 



n 


S3 


S4 


S30 


S40 


Sao/ Si,e 


-3 


34/7 


60712/1323 


26/7 


12088/441 


3022/1521 


-2 


27/7 


36457/1323 


19/7 


6019/441 


6019/3249 


-1 


20/7 


18376/1323 


12/7 


2008/441 


251/162 





13/7 


6469/1323 


5/7 


55/441 


11/45 


1 


6/7 


736/1323 


-2/7 


160/441 





Juszkiewicz et al. (1993) extend this analysis to the behavior of S30 for a Gaussian 
filter showing a good agreement with the theoretical predictions. I present here 
a new set of results for the skewness and the kurtosis for the density field using 
a top-hat filter to check, in particular, the accuracy of the scale dependence of 
the theoretical results for a CDM power spectrum. The measures have been made 
in an adaptive P 3 M simulation, kindly provided by Couchman (Couchman 1991), 
with CDM initial conditions, Ho = 50 km s _1 Mpc -1 and O = 1. It is done in a 
cubic box of side 400 Mpc with periodic boundary conditions and contains 2.1 10 6 
particles and at the end of the simulation the amplitude of the density fluctuations 
is cr 8 = 0.97. The resulting field has been filtered with a top-hat window function at 
four different radius and at two different timesteps corresponding to an expansion 
factor of 0.6 and 1 (in units of the final expansion factor). This procedure allows to 
check the stability of the result for different values of a. For each radius it is possible 
to calculate the index n of the power spectrum (Eq. [69]) as well as the parameter 
72 (= — dn/d log(-Ro)) and then the expected values of S3 and S4. These numbers 
are given in respectively the second to the fourth column of table 2. In parentheses 
is also given the value of S4 when the 72 term has been dropped. The remaining 
columns correspond to the measurement of these quantities together with a at the 
two different timesteps. 

The density has been measured in spherical cells of given radius centered on 
the 50 3 points of a grid. The actual second, third and fourth moments of the density 
distributions thus obtained have been calculated, as well as the ratios S3 and S4. 
The estimations of the errors have been made by taking the extreme deviations that 
are obtained when the volume of the simulation is divided in eight different parts. It 
provides an estimation of the sampling errors for the measurement of the skewness 
and the kurtosis expected in a cubic box of size 100/i _1 Mpc (h = H /100 km s _1 
Mpc -1 ) which is the order of magnitude of the size of available catalogues. As it can 
be seen the measurements are in good agreement with the theoretical predictions 
even for values of a slightly higher than unity. Unfortunately the determination of 
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the same quantities for the velocity field reveals, for a top-hat filter, to be more 
complicated and there are no available results. The results for the density field 
are however quite encouraging since if they are transposed to the divergence of 
the velocity field, it seems possible to measure the skewness and the kurtosis with 
enough accuracy to put strong constraint on O, for filtering radius up to 15 h~ x 
Mpc. At such a scale the perturbative calculations should be valid (and for the 
density field the analytical results are valid down to 5 /i _1 Mpc). 

Table 2: Comparison of the theoretical predictions with numerical 
simulations for a CDM initial spectrum, = 1 and Hq = 50 km s _1 
Mpc -1 . In parentheses is given the value of S4 when the 72 parameter is 
neglected, eq. (67). 





predictions 


measures a/a ~ 0.6 


a/ao = 1 


Ro/Mpc n 


72 £ 3 £4 


S3 S 4 


S3 S4 



10 - 1.34 - 0.46 

20 - 0.98 - 0.57 

30 - 0.73 - 0.64 

40 - 0.54 - 0.71 



3.2 17.7 (18.1) 

2.9 13.3 (13.7) 

2.6 10.6 (11.0) 

2.4 8.7 (9.2) 



0.92 3.2±g;| 18^ 

0.46 2.8jlg;J 13±ji 

0.29 2.4±J;£ 

0.21 2.2l} 5 16+1 



1.52 SAt i 20± 

0.74 2.8±° 5 6 13±" 

0.47 2.7±§;J 13± 

0.33 2.6±?;jj 15± 



A lot of numerical works still remain to be done to check the various 
features expected for the skewness and the kurtosis with regards to the cosmological 
parameters. Moreover it is not clear whether the validity domain is the same for 
any power spectrum and whether the errors associated with the finite size of the 
sample depend only on the local index of the power spectrum or on its global shape. 
In the latter case the sampling errors would have to be determined by numerical 
simulations for any model to be tested against the observations. 



4. COMMENTS 



4.1. Higher order moments 

Obviously the skewness and the kurtosis do not carry all the information 
available from the distribution of the divergence of the velocity field. The (connected 
part of the) moments of higher order also contain extra independent informations 
on the distribution function. The derivation of these moments requires similar 
perturbative calculations, and have been calculated (Bernardeau 1992a and b) 
when the filter function is neglected (which turns out to be equivalent to the limit 
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n — > —3) . The results are of the following form 

(S») c ~ 5 p (5 2 > p - 1 (71) 

for the density field and 

(<P) c ~ S^ey- 1 (72) 

for the divergence of the velocity field. The O and A of the coefficients S p are very 
weak but for the coefficients S p q they read 

S pe (A, O) = [/(A, fi)] 2 ^ S pe (A = 0, SI = 1) » qo.6(p-2) = 0, = 1). (73) 

This dependence with the cosmological parameters is expected to hold even with 
the application of a filter function (as verified for p = 3,4). It implies that whatever 
the number of moments you can measure it is impossible to give any constraint on 
A. Actually this property comes from the fact that the divergence of the velocity, 
at any order in the expansion (27), always depends on O and A like f(Q,A), 
9^(R ) oc f(Q, A), so that, at large scale, any statistical indicator can only give a 
constraint on this combination of O and A. 



4.2. The Zel'dovich approximation 



It has been claimed that the Zel'dovich approximation (Zel'dovich 1970) was 
an extremely good and powerful approximation for the large scale dynamics (e.g., 
Kofman et al. 1993). This approximation, based on a linear approximation for 
the displacement field, is a simplified dynamics that however reproduces most of 
the large-scale features of the gravitational dynamics. Obviously the results that 
have been obtained in the previous section do take into account corrections to 
the Zel'dovich approximation when necessary, that is the second or third order of 
the displacement field (see Bouchet et al. 1992 for an interesting discussion on 
this subject). It is anyway possible to calculate the large scale behaviors of the 
density and velocity moments in the simplified Zel'dovich dynamics. We obtain for 
P(k) oc k r 

4 -(3 + ?i); 

' [2-(3 + n)\; 



Zcl 



cZcl 

^38 



SP 1 

Q Zel 



272 50. . 7, n , 2 
- — (3 + n) + -(3 + n) 2 ; 



(74) 



9 
1 



n 



1.2 



f (3 + n) + 1(3 + n) 2 



The Q and A dependences are roughly the same as for the general results. As can be 
seen from these results the Zel'dovich approximation leads to an underestimation 
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of the correct values, especially for the velocity field since for instance for n = — 1 
the predicted skewness and kurtosis are zero. The use of a reconstruction method 
starting with the velocity field using the Zel'dovich approximation (Nusser & Dekel 
1993, Gramann 1993) to constrain O then may lead to a wrong estimation of O. 



4.3. The biases 

As mentioned previously the observations in galaxy fields provide good 
indications in favor of Gaussian initial conditions. However in case of biases 
between the galaxy distribution and the matter field, one can wonder whether the 
demonstration still holds. In fact, as pointed out by Fry & Gaztahaga (1992), the 
scaling between the moments of the galaxy distribution at large scale is expected 
to hold for a broad range of bias schemes. However, neither £3 or £4 can be used 
to give constraints on the galaxy linear bias, b (Eq. [3]). The reason is that the 
galaxy field skewness and kurtosis would be not only sensitive to the linear term b 6 
but also to the quadratic or cubic corrections. Following Fry & Gaztanaga, if we 
assume that, 

S 9 = bS+ b ^+ b ^ + ... (75) 



then we have, 



Q gal _ S 3 ^ 2 

^ = 1 + 12^ + 4^ + 4. 



(76) 



It shows that any measurement of high order moments involves a new parameter in 
the expansion (75), so that no new constraints on b alone can be given. However 

a departure between the measured values of and Sf 3 ^ and the theoretical 

predictions is the manifestation of biases that is 5 g is not equal to 5. For instance 
neither the skewness measured in the IRAS galaxy sample by (Bouchet et al. 
1993) nor the kurtosis are in agreement with the theoretical predictions. These 
are indications of the existence of biases in the galaxy distribution, although the 
value for b defined as the linear bias cannot be inferred from these measurements. 



4.4. A universal function to test the gravitational instability scenario 

Unlike the density field, the velocity divergence field is thought to be free of 
unknown biases. This supposes that it is possible to do a volume weighted filtering 
of the field as done by Juszkiewicz et al. (1993b) in a numerical simulation. For 
real data, however, it may not be easy (due for instance to undersampling in voids) 
but it is not the purpose of this paper to test the reliability of these tests on the 
present data. The skewness has already been proposed as a good indicator for the 
measurement of Q and measured in observational data set (Bernardeau et al. 1993). 
In fact it turns out that the kurtosis could be used as well to measure giving the 
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possibility to measure Q with two independent methods that are based on the same 
physical assumptions. Obviously the two results should be in agreement with each 
other. In fact there is a combination of the first three moments of the divergence 
of the velocity field that is almost independent of the cosmological parameters, 



The values taken by this ratio are given in table 1. The case n = 1 has not been 
given since the velocity skewness vanishes for n ~ 0.7 leading to an unstable ratio. 
In practice one should use the index observed in the velocity field rather than the 
one of the galaxy clustering since the latter may be affected by biases. However in 
the absence of such determination and to illustrate this result we can use the index 
observed at 10 h~ 1 Mpc in galaxy catalogues. It equals approximately -1 (there 
are numerous determinations but see for instance Peacock 1991) so that we expect 
a ratio of the order of 1.5 (table 1 and Figs. 5-8). Then, if the measured large 
scale velocity Hows fail to reproduce such a number, they cannot originate from 
gravitational instabilities with Gaussian initial conditions. 

4.5. Conclusion 

I gave here a series of results for the dependence of the skewness and the 
kurtosis of the cosmic fields at large scales with the cosmological parameters. These 
results have been exactly derived using perturbation theory. For the comparison 
with observational data however, numerous effects have to be been taken into 
account such as the redshift space distorsion for the density field (see Bouchet 
et al. 1992 for the skewness), the inhomogeneous Malmquist bias for the velocity 
field. Anyway the analysis of the statistical properties of the large-scale cosmic 
fields looks promising and it should be possible to constraint the density of the 
universe or to check the validity of the gravitational instability scenario. 
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12088 
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f^(n + 3) + f(n + 3) 2 
[f-(n + 3)] 2 



(77) 
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APPENDIX A: Geometrical properties of the top hat filter function 



This appendix is devoted to the derivation of general properties of the top 
hat window function, Wpjj(x). It is defined by 



W TH (x) = 1 if |x| < R ; 
Wpjj(x) = otherwise; 

for a scale Rq. The Fourier transform of this function is then given by 



(Al) 



W TR (k Rq) = ^J-p (sm(kRo) - kR cos(kR )) (A.2) 



where k is the norm of k. 



A.l. General properties 

Let me define two functions used in the main part: 



p( kl ,k 2 ) = l + ^ (A3) 

and 

Q(ki,k 2 ) = l-%^- i AA ) 

^1 "'2 

Then if one considers two wave vectors ki and k 2 on which one may have to 
integrate, a property of interest concerning the integrations over their angular parts, 
dOi and d0 2 , have already been given by Bernardeau (1993, Appendix B). It reads 



(47T) 



l -- [ dO!d0 2 Q^k^WTHdki+k,! R Q ) 

2 (A5) 



= 3 W TR (k 1 R Q )W TR (k 2 Ro). 



For the purpose of this paper we need a few other properties involving other 
geometrical dependences. 

They are the following: 



(47T) 



^2 ! dOidO, ^(ki.kajWTHdkx + kal Rq) 

= W TR (h Rq) W TR (k 2 R ) + -k 2 R W^ R (k 2 Rq) 
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and 



(4tt) ; 



J dO l( i0 2 Q(k 1; k 2 ) |ki + k 2 \R (|k! + k 2 | i? ) 



2 
3 



(fci # Wi- H (fci i2o) W XH (A; 2 i? ) + ^ TH (/ei i? ) fc 2 # W^ R (k 2 R ) 

(A.7) 



and 



(47T) 



\- J dOid0 2 P(k 1 ,k 2 )|k 1 +k 2 |i? ^ H (|k 1 + k 2 | i? ) 



=/ Cl R WTH'(h R ) W TR (k 2 R ) + -k x R W^ R (h R ) k 2 R W^ R (k 2 R ) 



~ 3%(^i Ro) k\ R\ W TR (k 2 R ). 

(A.8) 

The expression (A. 6) can be derived straightforwardly from the expression 
(A. 2) of the Fourier transform of the top hat window function. The integral over the 
angular part u = ki.k 2 //ci/c 2 can be done by the change of variable, x = |ki+k 2 | Rq, 
so that the integral reads, 



1 



(4tt) ; 



dOid0 2 



1 + 



h ki.k 2 

k 2 kik 2 

x dx 

syiki-kaiflo hk 2 R 2 



ki+k2| -Ro 



Wrndki + kal R ) 



i + 



k\k 2 



sinx — x cosx 



X"- 



(A9) 

The integration of (A. 9) can be done directly and leads to the expression (A. 6). 

The two other properties can be shown using the properties of the derivative 
of Wpjj(fc). Indeed we have 



1 



kWL R (k) = -W TR (k) + 



sin k 
~k~' 



(A10) 



The latter term can be written with a Bessel function, smk/k — \Ztt/2 Ji/ 2 (k) /k 1 / 2 , 
and and then decomposed as a sum over products of Bessel functions, 



sin 


ki + k 2 | 




k 


i +k 2 





m=0 



Jl/2+m(kl) Jl/ 2 + m (k 2 ) 

— r , m{— U), 



where 



u = 



k\ /2 

ki-k 2 

k x k 2 



(All) 



(A12) 
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and P m is the Legendre polynomial of order m. The general properties of the 
Legendre polynomials then lead to the relationships, 



1 f 1 

duP m (—u) = l if m = 0, otherwise; 

l 



2 



1 

- / dw uP m ( — u) = — 1/3 if m = 1, otherwise; (A. 13) 

2 i-i 

1 Z" 1 

- / dtt -u 2 Pm(— it) = 1/3 if m = 0, 2/15 if m = 2, otherwise. 

2 i-i 



Using these relationships we obtain, 
1 f 1 



du 



i , k± ■ 

1 + -z-U 

k 2 



sin ki + k 2 


2 sin k± sin k 2 


7T ^5/2 (&l) ^5/2(^2) 


|ki +k 2 | 


3 k\ k 2 


3 ^ 1/2 )fc 2 1/2 


sin ki + k 2 


sin fci sin k 2 


71" fcl ^3/2(^l) <^3/2(^2 


|ki +k 2 | 


fei /c 2 


2fc 2 fc ;/ 2 



It can also be noticed that , 



J 3/2 (k) = k 3 / 2 W TR (k)^ (A16), 
h/2{k) = -k 3 / 2 W^ R (k) (A17), 



so that 



sin 


ki +k 2 | 




k 


1 +k 2 





2 / sin k\ sin k 2 1 

3 V k\ k 2 



1 

2 









1 + T— U 




^2 



sin 


ki + k 2 


sin fci sin k 2 




k 


1 +k 2 


h k 2 



- ^iW^jjCfei) W TH (^)) (A18) 



-fc?W TH (fc 1 )Wr rH (* i 04.19) 



The relations (A. 7) and (A. 8) can then be deduced from (A.10), (A. 18) and (A.19). 
A. 2. Applications 

As an application of the previous results I give the value of few integrals of 
interest for the large scale cosmic fields. Let me define a 2 (Ro) by 

a2{Ro) = J J^ P{k) W ^ {k Ro) (A20) 
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where P(k) is the power spectrum. We can then notice that 

R o^^- = 2 J -^P(k)kR Q W TR (kRo)W^ R (kR ). 
The second derivative can also be calculated, 



Ro 



dR 



Ro 
+ 2 
+ 2 



da 2 (R ) 



dR 
■ d 3 k 
(27)" 3 
d 3 k 



2 / S^P(k) IkRoW^ikRo) 



(2tt) 3 

■P(fc) [kR ] 2 W TR (k R Q )W^ R (k Ro) 
P(k)kR W TR (kRo)W^ R (kR Q ). 



The expression (A. 22) can be simplified by the following property, 



so that 



Rq 



dR 



Ro 



k 2 W^ H (k) = -4kW^ R (k) - k 2 W TR (k), 
da 2 (R ) 

' =Z J (27)3 ^ [K) [ KTX ° V¥r i 

+ 2k 2 a 2 (R ) 



dR 



2 / -^P(k) \kR W^ R (kR ) 



+ 3Ro 



where 



k 2 = 1 



d 3 k 

a 2 (R )J (2^ 
I consider now the integrals, 
d 3 k! d 3 k 2 



dR 

k 2 P(k)W^ R (kR ) 



h(Ro) = 



(2tt) 3 (2tt) 3 



V(k u k 2 ) P(h) P(k 2 ) 



(A.21) 



(A.22) 



(A.23) 



(A.24) 



(A.25) 



x WT H (|ki+k 2 | Ro)W TR (k 1 Ro)W TR (k 2 Ro) 

h{Ro) = 1 00 Q(kl ' k2)p(/ci)p(/C2) 



x W XH (|ki + k 2 | Rq) W TR (ki R )W TR (k 2 Ro) 

(A.26) 

The relationships (A. 3) and (A. 4) give, 



h(R ) =a 4 (R ) 1 + - 



1 Ro da 2 (R ) 



6a 2 (R ) dR 



(A.27) 



I 2 (R ) = -<j*(R ) 
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The knowledge of integrals involving three wave vectors are also required. To 
ease the presentation I introduce simplified notations, 



W i = W TH (k i Ro), Wij = W TH (\k i + k j \ Rq), 
W ljk = W TR (\ki + kj + k fc | Ro) ■ 

'Pi, j = ^(kj, kj), Pij,k = P{ki + kj, kfc), Pijk = P{kii kj + kfc). 



and similar notations for Q(kj, kj). Then we have the results, 



J^f^^f P{h) P{h) P(h) Wl Wl2 W23 W " Vl ' 2 V2 > 3 



d 3 kx d 3 k 2 d 3 k 3 
(2tt) 3 (2tt) 3 (2tt) 3 



^(i? ) ( 1 + + l 7l 2 



P{h) P(k 2 ) P(k 3 ) W 1 W 12 W 23 W 3 V h2 V 3 ,2 



/ 
/ 



^^^s p[kl) p(k2) p(h) Wl Wl2 w ' 23 W3 Vv v " 



3 3 3 

W^^W^ p{kl) p{h) p{h) Wl Wl2 W2S W3 Vl > 2 Q2 > 3 

= \a\R Q ) (l + i 7l 

d 3 ki d 3 k 2 d 3 k 3 



6 

P(h) P(k 2 ) P(k 3 ) W 1 W 12 W 23 W 3 V 2 ,i Q 2)3 



(2tt) 3 (2tt) 3 (2tt) 3 

= 3' ' V 6 



2 a 6 (R ) (l + l-yi 



I W^§^§^ p{kl) p{h) p{h) Wl Wl2 W23 Ws Ql ' 2 Q2 > 3 = ^ 6(jRo) ' 

(A28) 

where 

_ dlog[q 2 (iZ )] d 2 log[a 2 (fZ )] 
71 dlo gj R ' 1 dlog 2 ^ ' 
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Another series of useful integrals are given by, 

P(fcl) P(fc 2 ) P(fc 3 ) W 2 W 3 W 123 Pl2,3 7>1,2 



' d 3 ki d 3 k 2 d 3 k 3 



(2tt) 3 (2tt) 3 (2tt) 3 



a 6 (i?o) ( l + i 7l + ^7i 2 



d 3 ki d 3 k 2 d 3 k 3 



(2tt) 3 (2tt) 3 (2tt; 



\ P(h) P(k 2 ) P(k 3 ) W x W 2 W 3 W 123 V 3 ,12 Pl,2 

3 3 3 

/ W^W^^ p{kl) p{h) p{h) Wl W2 w * Wl23 Vl2 > 3 Ql > 2 
/ W^§^§^ p{kl) p{k2) p{h) Wl W2 W3 Wl23 Vs > 12 Ql > 2 



P(h) P(k 2 ) P(k 3 ) W x W 2 W 3 W 123 Q 12)3 Vi,2 



d 3 kx d 3 k 2 d 3 k 3 
(2tt) 3 (2tt) 3 (2T) 3 

= V(i? )( 1+^71 



6 

/ P(/C2) P(/C3) Wl ^ ^ Wl23 Sl2 ' 3 Cl ' 2 = ^ 6(jRo) - 

(A30) 
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FIGURE CAPTIONS 



Fig. 1 : Variation of D 2 (Eq. [39]) and E 2 (Eq. [43]) as a function of O. Two 
cosmological hypothesis are considered: the solid line if for A = and the dashed 
line for a flat universe, O + A/3H 2 = 1. 

Fig. 2 : Variation of D 3 (Eq. [40]) and E 3 (Eq. [44]) as a function of O. The 
line symbols are the same as in Fig. 1. 

Fig. 3 : Large-scale behavior of the skewness of the density field (top) and of 
the divergence of the density field (bottom) as a function of O. The lines correspond 
to various index of the power spectrum index: from top to bottom n = —3, —2, 
— 1, 0, +1 in each panel and the solid lines are for A = and the dashed lines for 
Q + A/3H 2 = 1. 

Fig. 4 : The same as Fig. 3 but for the kurtosis of the density and velocity 
fields. In the latter case the n = 1 case is not presented since S40 is close to zero. 

Fig. 5 : The same as Fig. 3 but for the combination of the first three moments 
that is shown to have a very weak dependence with the cosmological parameters for 
the two fields. 

Fig. 6 : The moments as a function of the power spectrum index. The left 
panels are for the density field and the right for the divergence of the velocity field. 
The top panels give the skewness, the middle give the kurtosis and the bottom 
present a combination of these moments (Eq. [77]). The solid lines correspond to 
O = 1, A = 0, the long dashed line to O = 0. 3, A = and the short dashed lines to 
O = 0.3, A/3H 2 = 0.7. 

Fig. 7 : The moments as a function of scale in case of a CDM power 
spectrum. The top panels give the local index n (Eq. [70], thick lines) and 
72 (= — dn/d log(-Ro)) (thin lines) as a function of the smoothing radius. The 
functions S3 and S4 are calculated according to Eqs. [56-59] (thick lines) and thin 
lines correspond to the approximation (67-68) when the power spectrum is locally 
approximated by a power law of index n. The other symbols are the same as in Fig. 
6. The perturbation theory is thought to be valid for smoothing radius greater than 
~ 10/i -1 Mpc (h = -ffo/100 km s _1 Mpc _1 ). The squares with the error bars are 
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the measurements obtained in a CDM numerical simulation with = 1 presented 
in table 1 (and for a /do = 1). 

Fig. 8 : The moments as a function of scale in case of the power spectrum 
proposed by Peacock (1991) to fit the observed galaxy correlation function. Symbols 
are the same as in Figs. 6-7. 



